      subroutine solnem

      use param

#include <global.h>
      include 'ntbytes.h'
      include 'geom.h'
      include 'geomf.h'
      include 'xsec.h'
      include 'setls.h'
      include 'eigv.h'
      include 'nem.h'
      include 'sanm.h'

      real(NBF) ji11,ji12,ji21,ji22,jo11,jo12,jo21,jo22 &
     ,         fm1,fm2,sm1,sm2    

#define srcup
#ifdef srcup
      do im=1,ng
        do j=1,fnode
          nemsrc(j,im)=nemsrc(j,im)*rhzf(j)
          do k=1,4
            if (qsf(0,j,im).ne.0) then
!              qsf(k,j,im)=nemsrc(j,im)*(qsf(k,j,im)/qsf(0,j,im))
!              qsf(k,j,im)=0
            endif              
          enddo        
          qsf(0,j,im)=nemsrc(j,im)
        enddo
      enddo
#endif

!#define test
#ifdef test
      jinl(1,2,2)=0.29295
      jinl(2,2,2)=-0.04628
      jinr(1,2,2)=0.40886
      jinr(2,2,2)=-0.07616
      qs(0,2,2)=0.122425
      qs(1,2,2)=0.055205
      qs(2,2,2)=0.0118998
      qs(3,2,2)=0.0007475
      qs(4,2,2)=0.00460933
#endif
      do iin=1,3              
        do im=1,ng
          do i=1,fnode
            bt1=betaf(1,i,im)
            bt2=betaf(2,i,im)
            rxj=1/((25+384*bt2+96*bt1*(7+96*bt2))*(5+128*bt2+32*bt1*(7+160*bt2))) 
            rx0=1/(5+128*bt2+32*bt1*(7+160*bt2))
            rxf=1/(25+384*bt2+96*bt1*(7+96*bt2))
            f1=7+96*bt2
            f2=7+160*bt2
            f3=5+128*bt2
            f4=25+384*bt2
            f5=7+128*bt2
            f6=1+24*bt1
            f7=1+40*bt1
            f8=5+48*bt2
            f9=-5+3072*bt1*bt2
            f10=5+84*bt1
            f11=5+128*bt2
            f12=5+224*bt1
            f13=25+384*bt2
            f14=25+672*bt1
!
            t1=-rxj*(3072*f1*f2*bt1*bt1-f3*f4)
            t3=-768*rxj*bt1*(24*f5*bt1+96*bt2+5)
            t5=-256*rxj*bt2*(24*f5*bt1+96*bt2+5)
            t7=rxj*(150528*bt1*bt1+8960*bt1-49152*f6*f7*bt2*bt2+125)
!
            t2=-64*rxj*bt1*(256*f8*bt2+35)
            t4=192*f9*rxj*bt1
            t6=64*f9*rxj*bt2
            t8=-256*rxj*(64*f10*bt1+5)*bt2
!
            t9=20*f3*rx0*bt1
            t10=-640*rx0*bt1*bt2
            t11=-1920*rx0*bt1*bt2
            t12=20*f12*rx0*bt2
!
            t13=30*f4*rxf*bt1
            t14=-2880*rxf*bt1*bt2
            t15=-8640*rxf*bt1*bt2
            t16=30*f14*rxf*bt2
!
            a011=2*t9*rhzf(i)+xsrf(i,im)
            a012=2*t11*rhzf(i)-2*xsrf(i,im)
            a021=2*t10*rhzf(i)-2*xsrf(i,im)*de3
            a022=2*t12*rhzf(i)+de3*(4*xsrf(i,im)+5*xstf(i,im))
!
            as11=-7*t9*rhzf(i)
            as12=-7*t11*rhzf(i)
            as21=-7*t10*rhzf(i) 
            as22=-7*t12*rhzf(i) 
!
            b011=-(t1+t2-1)*rhzf(i)
            b012=-(t3+t4)*rhzf(i) 
            b021=-(t5+t6)*rhzf(i)
            b022=-(t7+t8-1)*rhzf(i)
!
            af11=(-de5*(56*t13+24*t14)+15)*sigdf(1,i,im)+xsrf(i,im)
            af12=-de5*(56*t15+24*t16)*sigdf(1,i,im)-2*xsrf(i,im)
            af21=-de5*(8*t13+32*t14)*sigdf(2,i,im)-2*de3*xsrf(i,im)
            af22=(-de5*(8*t15+32*t16)+15)*sigdf(2,i,im)+de3*(4*xsrf(i,im)+5*xstf(i,im))
!
            bf11=4*de5*(7*t1-7*t2+3*t5-3*t6+7)*sigdf(1,i,im)
            bf12=4*de5*(7*t3-7*t4+3*t7-3*t8+3)*sigdf(1,i,im)
            bf21=4*de5*(t1-t2+4*t5-4*t6+1)*sigdf(2,i,im)
            bf22=4*de5*(t3-t4+4*t7-4*t8+4)*sigdf(2,i,im)
!
            tt1=af22*bf11-af12*bf21
            tt2=af22*bf12-af12*bf22
            tt3=af21*bf11-af11*bf21
            tt4=af21*bf12-af11*bf22
            rr1=1/(af11*af22-af12*af21)
            rr2=1/(af12*af21-af11*af22)
            qt1=de3*qsf(1,i,im)+de5*qsf(3,i,im)
            qt2=de5*qsf(2,i,im)-3*de35*qsf(4,i,im)
            ji11=jinl(1,i,im)
            ji12=jinr(1,i,im)
            ji21=jinl(2,i,im)
            ji22=jinr(2,i,im)          
! First moments
            fm(1,i,im)=rr1*(tt1*(ji12-ji11)+tt2*(ji22-ji21)+de3*qt1*(2*af12+3*af22))
            fm(2,i,im)=rr2*(tt3*(ji12-ji11)+tt4*(ji22-ji21)+de3*qt1*(2*af11+3*af21))
            fm1=fm(1,i,im)
            fm2=fm(2,i,im)
!
            s011=sigdf(1,i,im)*(de25*(392*t9+168*t10)-7)
            s012=sigdf(1,i,im)*de25*(392*t11+168*t12)
            s021=sigdf(2,i,im)*de25*(56*t9+224*t10)
            s022=sigdf(2,i,im)*(de25*(56*t11+224*t12)-7)
!
            ss11=sigdf(1,i,im)*(-de25*(1372*t9+588*t10)+35)+xsrf(i,im)
            ss12=-sigdf(1,i,im)*de25*(1372*t11+588*t12)-2*xsrf(i,im)
            ss21=-sigdf(2,i,im)*de25*(196*t9+784*t10)-2*de3*xsrf(i,im)
            ss22=sigdf(2,i,im)*(-de25*(196*t11+784*t12)+35)+de3*(4*xsrf(i,im)+5*xstf(i,im))
!
            cs1=qsf(0,i,im)-rhzf(i)*((-1+t1+t2)*(ji11+ji12)+(t3+t4)*(ji21+ji22))
            cs2=-2*de3*qsf(0,i,im)-rhzf(i)*((t5+t6)*(ji11+ji12)+(-1+t7+t8)*(ji21+ji22))
            cs3=0.2*qsf(2,i,im)-3*de35*qsf(4,i,im)-1.12*sigdf(1,i,im) &
               *((7*(1+t1+t2)+3*t5+3*t6)*(ji11+ji12)+(3*(1+t7+t8)+7*t3+7*t4)*(ji21+ji22)) 
            cs4=-2*de15*qsf(2,i,im)+2*de35*qsf(4,i,im)-1.12*sigdf(2,i,im) &
               *((1+t1+t2+4*t5+4*t6)*(ji11+ji12)+(4+t3+t4+4*t7+4*t8)*(ji21+ji22))
!
            ta=s012*s021-s011*s022
            tb=s011*ss21-s021*ss11
            tc=s022*ss11-s012*ss21
            td=s021*ss12-s011*ss22
            te=s012*ss22-s022*ss12
            tf=ss12*ss21-ss11*ss22
            tg=a022*ss21-as21*s022
            th=as22*s022-a022*ss22
            ti=as21*ss22-as22*ss21
            tj=as21*s012-a022*ss11
            tk=a022*ss12-as22*s012
            tl=as22*ss11-as21*ss12
            tm=a021*ss21-as21*s021
            tn=as22*s021-a021*ss22
            tp=as21*ss22-as22*ss21
            tq=as21*s011-a021*ss11
            tr=a021*ss12-as22*s011
            ts=a021*as12-a011*as22
            tt=a011*ss22-as12*s021
            tu=as22*s021-a021*ss22
            tv=a011*ss12-as12*s011
            tw=a021*as11-a011*as21
            tx=a011*ss21-as11*s021
            ty=a011*ss11-as11*s011
            tta=a011*a022-a012*a021
            ttb=a012*s021-a011*s022
            ttc=a021*s022-a022*s021
            ttd=a012*s011-a011*s012
            tte=a021*s012-a022*s011  
! Even moments
            denom1=1/(as12*(as21*ta+a022*tb+a021*tc)+as11*(-as22*ta+a022*td+a021*te) &
                  +a012*(-as22*tb-as21*td+a021*tf)-a011*(as22*tc+as21*te+a022*tf)) 
           denom2=1/(-(ts*s012+a022*tv-a012*tr)*(tw*s022+a022*tx-a012*tm) &
                  +(tw*s012+a022*ty+a012*tq)*(ts*s022+a022*tt+a012*tn))
           pphi(1,i,im)=(-cs1*(as22*tc+as21*te+a022*tf)+cs2*(as12*tc+as11*te+a012*tf) &
                         +cs3*(as12*tg+as11*th+a012*ti)+cs4*(as12*tj+as11*tk+a012*tl))*denom1          
           aflx(2,i,im)=-(cs1*(as22*tb+as21*td-a021*tf)+cs2*(-as12*tb-as11*td+a011*tf) &
                          +cs3*(as12*tm+as11*tn+a011*tp)+cs4*(as12*tq+as11*tr+a011*tl))*denom1 
           aflx(1,i,im)=pphi(1,i,im)-2*aflx(2,i,im)
           sm(1,i,im)=(cs1*(-as22*ta+a022*td+a021*te)+cs2*(as12*ta-a012*td-a011*te) &
                      +cs3*(s022*ts+a022*tt+a012*tu)-cs4*(s012*ts+a022*tv-a012*tr))*denom1   
           sm(2,i,im)=((cs1*ttc+cs2*ttb+cs4*tta)*(s012*tw+a022*ty+a012*tq) &
                      -(cs1*tte+cs2*ttd+cs3*tta)*(s022*tw+a022*tx-a012*tm))*denom2 
            sm1=sm(1,i,im)
            sm2=sm(2,i,im)
! Outgoing partial currents
            joutl(1,i,im)=t1*ji11+t3*ji21+t2*ji12+t4*ji22+t9*pphi(1,i,im)+t11*aflx(2,i,im) &   
                         -t13*fm1-t15*fm2-3.5*sm1
            joutl(2,i,im)=t5*ji11+t7*ji21+t6*ji12+t8*ji22+t10*pphi(1,i,im)+t12*aflx(2,i,im) &    
                         -t14*fm1-t16*fm2-3.5*sm2
            joutr(1,i,im)=t1*ji12+t3*ji22+t2*ji11+t4*ji21+t9*pphi(1,i,im)+t11*aflx(2,i,im) &    
                         +t13*fm1+t15*fm2-3.5*sm1
            joutr(2,i,im)=t5*ji12+t7*ji22+t6*ji11+t8*ji21+t10*pphi(1,i,im)+t12*aflx(2,i,im) &    
                         +t14*fm1+t16*fm2-3.5*sm2
            jo11=joutl(1,i,im)
            jo21=joutl(2,i,im)
            jo12=joutr(1,i,im)
            jo22=joutr(2,i,im)
! Surface flux
            sflx(1,i,im)=0.32*(7*(ji11+jo11)+3*(ji21+jo21))
            sflx(2,i,im)=0.32*(ji11+jo11+4*(ji21+jo21))
            if (i.eq.fnode) then
              sflx(1,i+1,im)=0.32*(7*(ji12+jo12)+3*(ji22+jo22))
              sflx(2,i+1,im)=0.32*(ji12+jo12+4*(ji22+jo22))
            endif    
            sfl(1,i,im)=0.32*(7*(ji11+jo11)+3*(ji21+jo21))
            sfr(1,i,im)=0.32*(7*(ji12+jo12)+3*(ji22+jo22))
            sfl(2,i,im)=0.32*(ji11+jo11+4*(ji21+jo21))
            sfr(2,i,im)=0.32*(ji12+jo12+4*(ji22+jo22))

            call updjin(i,im)  

          enddo  ! for i
        enddo  ! for im

        call updsrc
! Net current
        do im=1,ng
          do i=1,fnode
            jnetf(1,i,im)=jinl(1,i,im)-joutl(1,i,im)
            jnetf(2,i,im)=jinl(2,i,im)-joutl(2,i,im)
            pphi(2,i,im)=aflx(2,i,im)
          enddo
          jnetf(1,i,im)=joutr(1,fnode,im)-jinr(1,fnode,im)
          jnetf(2,i,im)=joutr(2,fnode,im)-jinr(2,fnode,im)
        enddo
      enddo  ! for iin 

      call dumflx
            
      end subroutine

      subroutine updsrc

      use param

#include <global.h>
      include 'ntbytes.h'
      include 'geom.h'
      include 'geomf.h'
      include 'xsec.h'
      include 'setls.h'
      include 'eigv.h'
      include 'sanm.h'
      include 'nem.h'      

      real(NBF) fmt1,fmt2,smt1,smt2

! 5 coefficients of the 4-th order flux
      do i=1,fnode
        fscf=0
        do im=1,ng
          fmt1=fm(1,i,im)
          fmt2=fm(2,i,im)
          smt1=sm(1,i,im)
          smt2=sm(2,i,im)

          cfp10=pphi(1,i,im)
          cfp11=0.5*(sfr(1,i,im)-sfl(1,i,im))
          cfp12=pphi(1,i,im)-0.5*(sfl(1,i,im)+sfr(1,i,im))
          cfp13=5*(fmt1-de3*cfp11)
          cfp14=7*de3*(cfp12-5*smt1)
          cfp20=aflx(2,i,im)
          cfp21=0.5*(sfr(2,i,im)-sfl(2,i,im))
          cfp22=aflx(2,i,im)-0.5*(sfl(2,i,im)+sfr(2,i,im))
          cfp23=5*(fmt2-de3*cfp21)
          cfp24=7*de3*(cfp22-5*smt2)
!
          flp(0,im)=cfp10-2*cfp20
          flp(1,im)=cfp11-2*cfp21
          flp(2,im)=cfp12-2*cfp22
          flp(3,im)=cfp13-2*cfp23
          flp(4,im)=cfp14-2*cfp24

! Update fission source coefficients
          do l=0,4
            fscf(l,i)=fscf(l,i)+xsff(i,im)*flp(l,im)
          enddo
        enddo
! Update source coefficients
        do im=1,ng
          sscf=0  
          do m=1,ng   
            do l=0,4
              sscf(l,i)=sscf(l,i)+xsmf(i,im)%from(m)*flp(l,m)      
            enddo
          enddo  
          do l=0,4
            qsf(l,i,im)=reigv*xchif(i,im)*fscf(l,i)+sscf(l,i)          
          enddo
        enddo
      enddo

      return
      end subroutine
